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The transport of material by spatio-temporally complex flow fields has widespread 
application to geophysical and industrial processes. In recent years, advances 
in methods to detect Lagrangian Coherent Structures (LCS), these being key 
material lines within the flow field that organize flow transport, have enabled 
significant breakthroughs and profound new insight. One of the earliest devel¬ 
oped LCS detection methods employs the finite-time Lyapunov exponent held 
(FTLE), and while newer and more rigorous methods exist, the simplicity of 
the FTLE approach, both in its implementation and its identification of active 
regions of a flow field, means that it is still widely used. As such, the thrust of 
this work is to highlight good practices in calculating the FTLE held, to provide 
a practical means to identify the most active material lines of a flow field that lie 
along maximal FTLE ridges, and to classify the types of deformation associated 
with an FTLE ridge. The methods we refine and advocate, which build on 
previous established results, are tested on an analytical model and a numerical 
model data set of an ocean surface flow. 


I. INTRODUCTION 

There has been substantial recent effort to develop, establish and advance rigorous, ob¬ 
jective methodologies and associated data analysis tools for understanding how advective 
transport is organized in complex, time dependent fluid flow j^ * ^^ *^. One of the first tools 
developed for this purpose was the Finite Time Lyapunov Exponent (FTLE), which reveals 
the regions of the flow held that undergo the greatest stretching for the time window con- 
siderecP^. More sophisticated and rigorous approaches have subsequently been developed; 
in particular, geodesic methods identify structures such as strainlines and shearlines, whose 
properties are very well defined and for which there are numerically stable tool^. The 
progress has been substantial, particularly for two-dimensional flows that are relevant to 
such important scenarios as the evolution of an oil spill on the ocean surfac#^. 

Although more advanced methods now exist, obtaining the FTLE held and visualizing 
the FTLE ridges remains a popular and insightful means for investigating the organization of 
transport in complex flows. A healthy degree of caution and understanding are needed when 
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interpreting the results. For example, it should be recognized that since FTLE ridges are 
simply representations of the time history of what happens to material elements, it may be 
that an entire FTLE ridge exists only because the material elements of which it is comprised 
are transported past an important local Lagrangian flow feature (e.g. a hyperbolic core)?®. 
It is unarguably the case, however, that FTLE ridges represent the most kinematically 
active (i.e. most local stretching) material lines of the flow held, and this may be the most 
important consideration for a study. Furthermore, strainlines and shearlines, these being 
key material lines identihed by the geodesic approach^, have no obligation to align with 
FTLE ridges in general, compressible two-dimensional hows, and so one cannot assume that 
an FTLE ridge will be so marked (although it is always worth checking). 

While FTLE analysis has been, and continues to be, widely utilizecP^®, there is scope 
for improvement. For example, as with any numerical method, it is important to check 
convergence of the FTLE values with system parameters, such as the velocity held resolution 
or the cluster size used for hnite-diherence approximation; while somewhat rudimentary, 
this topic has yet to be systematically addressed in the FTLE literature. It also remains a 
practical challenge to robustly extract ridges of the FTLE held, which provides a simpler 
presentation of the results and enables material line advection, thereby revealing how ETLE 
ridges evolve and shape transport. Finally, an outstanding issue is that FLTE ridges do not 
reveal what type of local deformation they represent (i.e. normally hyperbolic repulsion, 
Lagrangian shear or tangential stretching) and there has been only one previous stud}^ 
that attempted to take this further step; the ability to do so would provide more insight 
into the nature and signihcance of an FTLE ridge. 

In this paper, recognizing the continued use by many of FTLE analysis particularly when 
studying geophysical hows, we revisit the topic and seek to address the aforementioned out¬ 
standing issues. In section [IT], we demonstrate good practices in checking convergence of the 
FTLE held, as well as outlining an alternative method for calculating FTLE helds that is 


highly accurate for analytic data sets. Next, section III presents a practical extraction and 
rehnement scheme for determining FTLE ridges with sufficient accuracy that, despite being 
the most sensitive features in the how held, they can be faithfully advected. Section [TV] then 
attempts to classify FTLE ridges according the type of local deformations they represent, 
identifying the inherent challenges and resolving what classihcations are reasonably achiev¬ 
able and under what circumstances. Having addressed these issues using an analytic model 
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as a test case, in section |V] we proceed to apply our findings to a real-world case study 
of an ocean surface flow, using a data set generated by a high quality numerical model. 


Finally, in section VI, we present our conclusions. Some technical details are presented in 
the appendices. 


II. FTLE CALCULATIONS 

The principal task underlying the determination of the FTLE held is the calculation 
of the right Cauchy-Green strain tensor (henceforth referred to as the CG tensor) for the 
physical domain and time window of interest. In this section, we discuss the standard hnite- 
difference methods and a less well utilized approach to calculate the GG tensor, investigate 
the impacts of discretization and noise, and discuss good practices. While the focus of the 
present study is FTLE ridges, it should be noted that the following calculations are also 
necessary for strainline and shearline based methods. 


A. Methods 

The hrst step in calculating the GG tensor is to calculate the flow map, which maps a fluid 
element from its initial position Xq = (01,02) at time to to its final position x = {xi,X2) 
at time t, and is represented as = x{t;xo,to). To elucidate the notation, the 

initial condition, also the Lagrangian coordinates, (01,02) are used to differentiate from the 
Eularian coordinates {xi,X2). Given a velocity field u{x,t), solutions to the flow map are 
determined by solving the system of Ordinary Differential Equations (ODEs): 

dx . . .. 

w 

for an initial position of Xq at time to- We use the Matlab solver ode 45 , which is a variable 
step Runge-Kutta 4 *^ and 5 *^ order ODE solver, for all calculations in this paper. If the 
velocity field is analytic, then velocities can be calculated directly. If it is a discrete velocity 
held data set, however, velocities at off-grid locations required by the ODE solver must be 
obtained by interpolation, the impacts of which are investigated later in this section. 

In order to map a set of initial points ..., the equations of motion ([^ 

must be solved for each initial condition. This process is accelerated by simultaneously 
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solving the series of equations: 




1 

x^^\x^'^\t) = u{x^^\t). 


( 2 ) 


Simultaneous advection takes advantage of the built in vectorization of Matlab and imposes 
the same variable time step sequence and numerical tolerances to all trajectories. For sen¬ 
sitive systems, failure to advect trajectories simultaneously may result in erroneous results, 
particularly if the relative tolerances are large. As good practice, the absolute and rela¬ 
tive tolerance for the solver should be systematically reduced so that convergence of the 
trajectories is achieved at the associated but necessary cost of increased computational run 
time. 

For a given initial condition Xq, the advected hnal position of a nearby point, Xq -|-e, can 
be approximated by: 

Fl(x„ + £) = Fl(x„) + VF‘Jx„)e + 0 (|ep), ( 3 ) 


where the second term in this expansion contains the flow map gradient: 


VFl{xo) = dxi/dttj 


for i,j 

Xq 


1 , 2 , 


( 4 ) 


Finite-difference methods have been the primary means for calculating the flow map gradi¬ 
ent, with nearest neighbors on a regular grid initially being used for the derivative calcu- 
lationd^. Uniform-grid based methods pose an issue when trying to improve the accuracy 
of the hnite-difference scheme, however, due to the rapidly increasing computational de¬ 
mand as the grid spacing is reduced. Unstructured mesh methods have been utilized, with 
dynamic mesh rehnement used to improve the resolution of flow map calculations in high 
FTLE region^l^l. While this improves accuracy, the need for mesh rehnement complicates 
the calculation, increasing the computational demand. The use of cluster^, which adds 
four off grid points of adjustable distance to calculate the hnite-difference term for each 
on-grid location, has become widely used. This method allows for Cartesian grids with run 
times independent of the desired accuracy, at the not unreasonable expense of a hve-fold 
increase in computational cost. Arbitrarily hne accuracy cannot be achieved as the cluster 
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size is reduced, of course, as the difference in advected position eventually becomes the same 
order as the numerical resolution, resulting in errors. Underlying all these hnite-difference 
techniques is the possibility that neighboring points may not remain nearby throughout 
advection. In such a situation, the hnite-difference approximation of a local deformation 
may itself not be valid. To compensate, one can renormalize the cluster around the central 
point throughout advection preventing nearby trajectories from becoming non-locaP^ this, 
of course, increases computational complexity and demands. 

An alternative, but relatively untested, approach for calculating the terms of the how 
map gradient is to use the advected-gradient method, which solves a system of ODEs that 
yields the terms of both K (x) and 'VF}^{xq) directly. This is achieved by solving Q and 
the following system of equations, obtained via the chain rule, as follows: 


d dxi d dxi dxk dui 


for i, j = 1 , 2 , 


( 5 ) 


dt dttj dttj dt dttj dxk 
where we have used standard index notation convention for summing repeated indices; the 
initial conditions for these equations are: 


xito) = Xo, 


dxi dx2 


= 1 , 


dxi dx2 


= 0 . 


( 6 ) 


dai da2 ’ da2 dai 

An advantage of the advected-gradient approach is that for an analytic system, where the 
velocity held and its gradients are completely known, this method will be accurate. The six 
coupled equations typically take longer to solve than the ten equations needed for a cluster 
based hnite-diherence approach (two equations for each of the hve particles in a cluster); 
this is due to the large and rapidly changing values of the how-map-gradient terms requiring 
more computation time to satisfy equivalent numerical tolerances. It also remains to be 
determined to what extent the advected-gradient method is compromised by discretization 
and noise. 

Once the how map gradient held has been determined, the FTLE held is readily extracted 
from the CG tensor held: 


C = IVF'l-IVF'l, 


( 7 ) 


where * corresponds to the transpose operator. By dehnition, the FTLE is: 

where A2 is the largest eigenvalue of C. 
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B. Results 


To test the FTLE methods, we use an autonomous analytic model that has analytic 
solutions for all the quantities we are concerned with, so that in all cases there is a ‘true’ 
result against which to compare the results of our various numerical calculations. Here we 
outline the model, more detail on which is given in appendix The autonomous analytical 
system is based on the nonlinear system of equations: 

X2=X2'-X2, (9) 

within the domain D = {Xi G [—1,1],X2 G [—1,1]}. To add spatial complexity to the flow 
held we introduce a coordinate transformation: 


Xi = Xi cos(r) — X 2 sin(r), 

=-^2 cos(r) + Xi sin(r), (10) 


where r = \/xi^ + X 2 ^ = a/Xi^ + X 2 ^. The equations of motion for the transformed coordi¬ 


nates are calculated by taking the time derivative of equation (10): 


xi = — xf — { - r){xl — xl)x^ cos(2r) 

-h ^ 0:2 - ^X2{3xI + xl) “ “ r)2xixl^ sin(2r), 

X 2 = (^X 2 -f- a ;2 + (- r)[x\ — x\)x^ cos(2r) 

-I- — ^Xi{x\ + 3 X 2 ) — — r)2xlx‘^ sin(2r). (11) 

Due to the transformation, the repelling invariant manifolds, X 2 = ±1, now correspond to 
the left and right boundaries. The origin remains a stationary point in phase space and 
is still a hyperbolic cor^^. Neither the original system nor the transformed system satisfy 
continuity, which is intentionally so in order to permit a greater variety of local deformations 
than possible for an incompressible system. While the level of compressibility in the model 
system is strong, in the ocean, for example, upwellings and downwellings can cause the 
surface velocity held (i.e. not accounting for the vertical velocity) to be highly compressible. 

The analytic FTLE held for a time window t — to = 2 is presented in hgure[^a). The 
largest FTLE values are along the left and right boundaries and through the center of 
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FIG. 1. (a) Analytic FTLE field for the autonomous system with t — to = 2. (b) Difference 
between the analytic and finite-difference FTLE fields, (c) Difference between the analytic and 
advected-gradient FTLE fields. 


the domain. While the largest values of the FTLE held exceed 2, there is a large portion 
of the domain where the FTLE value is negative, highlighting the fact that the system 
does not satisfy continuity. For comparison, the error in the FTLE helds calculated via 
the hnite-difference and advected-gradient approaches are presented in hguresj^b) and (c), 
respectively. The hnite-difference calculation was most accurate when using a cluster size 
6a = 10“® and a relative tolerance of 10“^; these parameters yielded a spatially averaged 
error on the order of 10“®, with the largest errors occurring along the curvilinear left and 
right boundaries. While the errors in the FTLE values along the boundary are both positive 
and negative, there is a persistent underestimation of the FTLE held across the center of the 
domain. The advected-gradient result has a similar structure to the hnite-diherence error 
held; the former is on average two orders of magnitude smaller, however, with a spatially- 
average error of order 10“^°. For both calculations, identical relative tolerances were used, 
therefore the additional error is solely due to the hnite-diherence approximation. 

Typically, for FTLE calculations the analytic velocity held is not known and instead a 
velocity held data set discretized in space and time is available. It has been established that 
FTLE ridges are robust to isolated errors in the velocity held either spatially or temporally 
given that the error is smalP, but persistent errors resulting from discretization or noise may 
not satisfy this condition. When analyzing ocean models, spatial resolution is fundamental 
to determining the scale of FTLE structures that are reliably detected. A number of stud¬ 
ies have investigated the ehects that ocean model resolution has on the resulting FTLE or 
hnite-size Lyapunov exponent held^^*^^*^. One particular study systematically subsampled 
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turbulent ocean models to study the impact of spatial and temporal on the FTLE calcula- 
tionCSl. Another study has investigated the impact of noise and discritization on an FTLE 
held calculated via hnite-difference on a uniform helci^, here we focus on the relative im¬ 
pacts of both on the different FTLE calculation methods. Our goal is to study the impact 
of varying spatial resolution and the presence of noise on the hnite-difference based FTLE 
calculation relative to the analytically calculated values. 

To investigate the ehects of discretization, we considered a uniform grid of data from 
our analytical model for the range of grid spacings Ax = 2“"^ to 2“^^, utilizing Matlab’s 
griddedinterpolants to accelerate the interpolation. Because Matlab’s interpolant is faster 
when performing simultaneous as opposed to sequential interpolations, the simultaneous 
advection scheme given by equation (|^ again has signihcant advantages over the advection 
of individual trajectories. To quantify the degree of error in regions of interest, we calculated 
the spatially averaged error of the FTLE held, <he, in regions where $ > 1. This measure 
is presented as a function of the velocity held resolution for both the hnite-diherence and 
advected-gradient methods in hgurej^a); furthermore, for the hnite-diherence method the 
results are obtained for cluster sizes 5a = 10“^, 10“®, 10“®, and 10“^°. Both methods show 
order 0.1% errors for the coarsest velocity held resolution. As Ax is reduced, the hnite- 
diherence results converge to the numerical errors obtained when working with the analytic 
velocity held. The error of the advected-gradient method, however, is typically two orders 
of magnitude larger than those of the hnite-diherence method. This is attributed to the 
sensitivity of the how map gradient terms to errors in the velocity gradient held, and reveals 
that the hnite-diherence method performs better for discretized data. 

In addition to discrete data sets, another concern is that velocity held data sets may 
be noisy. To systematically study the ehect of noise, a data set with spatial resolution 
Ax = was used (so that noise was the dominant factor over discretization) and at each 
grid point uniformly distributed noise of magnitude ||e|| was added to the velocity held, 
which is of order 10“^. The results for <Fe for the two methods are presented in hgurej^b). 
For large values of ||e||, the advected-gradient and hnite-diherence methods are similarly 
impacted by noise, with the former method somewhat outperforming the latter as noise 
levels are reduced. 

Based on these studies, the introduction of a discretized held, with or without noise, 
ehectively eradicates the greater accuracy of the advected-gradient method, and thus the 
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FIG. 2. (a) Spatially averaged FTLE error for regions where > 1 for the advected-gradient 
(dashed) and finite-difference methods (solid) calculated from a discrete velocity field with varying 
resolution, Ax; the cluster spacing, 5a, is also varied. All calculations used a relative tolerance 
of 10“^ for the ODE solver, (b) Error calculated from a discrete velocity field with a resolution 
Ax = with uniformly distributed noise of maximum magnitude ||e||. 

cluster-based finite-difference method should be employed. While this study has been per¬ 
formed for a simple analytic system, we expect these results to carry over to more spatially- 
complex, time-dependent discretized velocity helds, as we performed similar studies for the 
double-gyre flow that generated similar results. 


III. FTLE RIDGES 

The standout features of an FTLE held are the maximal ridges, which correspond to 
regions that undergo maximum Lagrangian separation of neighboring huid elements. As has 
been pointed out, FTLE ridges do not necessarily coincide with variationaP and geodesicP 
dehnitions of LCS, but they may nevertheless be considered a type of LCS, in the sense that 
they identify the most kinematically active material lines in a how held for a given time 
window. It is desirable to be able to accurately identify the material lines that reside along 
FTLE ridges in order to form a simplihed picture of the how transport, to enable these 
structures to be advected over the time window, and, as we shall see, to enable classihcation 
of the associated types of deformation. In so doing, it should be appreciated that being in 
such active regions of the how held, these material lines will likely be greatly deformed as 
they are advected and may not represent a barrier to transport, per se. 
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A. Methods 


To develop the foundation for detecting and classifying FTLE ridges, it is first necessary 
to define a ridge. For two-dimensional systems, ridges of the FTLE field are one-dimensional 
lines that mark generalized local maxima of the held. A “height ridge” is a line made up of a 
set of points that are locally maximum in the direction of the greatest curvature of the FTLE 
held projection®. While this dehnition lends itself to a point by point evaluation, it does not 
necessarily create connected lines forming ridges, and attempts to amend the height ridge 
dehnitioEpSI resulted in an over-constrained systempS. An alternative to the height ridge, 
the “watershed ridge” divides the system into disjoint regions based on global stationary 
points of the helcP. The ridges are identihed as slope lines, these being trajectories of the 
FTLE gradient held that connect the saddle points in the system to the local maxima. One 
variation of the watershed ridge removed the requirement that it start at stationary points 
of the FTLE held but imposes a requirement that the slope line be a normally attracting 
invariant manifolcP^l. Both the height ridge and watershed ridges dehnitions are practically 
challenging to implement in the presence of noisy data. In the case of the height ridge, the 
challenge arises when trying to calculate the Hessian matrix of the FTLE held. And while 
the watershed dehnition does naturally produce connected lines, noise may produce many 
stationary points in the FTLE held and signihcantly contaminate the gradient held. 


Building on the dehnitions of height ridges and watershed ridges, we form a numerically 
tractable ridge dehnition. Principally, a ridge should be a set of connected points that are 
generalized local maxima of the FTLE held, so points on the ridge will be a local maximum 
FTLE value in the direction normal to the ridge. Compared to the height ridge dehnition, 
this dehnition loosens the requirement that the ridge normal be aligned with the smallest 
eigenvector of the Hessian of the FTLE held, thus removing the need to calculate second 
derivatives of the numerically evaluated ETLE held. To form an actual line, ridges will be 
everywhere tangent to the gradient of the FTLE held, like the slope lines that form watershed 
ridges, but there is no requirement that the starting points of the ridges be stationary points 
of the FTLE held, which removes the need to accurately identify the numerically evaluated 
stationary points. With these criteria in hand, we dehne a normal-maximum ridge as the 
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line 7 where 'ix G 7 : 

V$ 

e(a?) = ( 12 a) 

n{x)*Hq,{x)n{x) < 0 , ( 12 b) 

where e{x) and n{x) are the tangent and normal to 7 at x, respectively. The hrst condition 
ensures that the ridges are slope lines, and the second condition ensures points are a local 
maximum in the direction normal to the ridge. In practice, the Hessian is not actually 
calculated, and instead FTLE values of points nearby and along a normal to the ridge are 
used to reasonably check the local maximum criteria. 

There already exists a numerical algorithm that with minor modihcations can be used 
to locate normal-maximum ridge^l^^. The basic ridge tracking algorithm, as illustrated in 
hgure is as follows: 

1. Locate seed points, x{s = 0), which are local FTLE maxima in the initial step direc¬ 
tion, V<h. To accelerate the search, points are selected from a grid, Q, made up of 
vertical and horizontal lines that divide the system. Start two trajectories from each 
seed point with initial step directions opposite of each other 

2. Step in this direction a distance As, and calculate the FTLE values at the new position 
as well as two other points, one either side in a direction normal to the step direction 
taken. 

3. Fit a parabola to the FTLE values calculated and add the corresponding location of 
the maximum, a;(s -|- As) to the ridge. 

4. Use the previous step position, a;(s), and the current position to calculate the tangent 
vector. 

5. Repeat steps 2-4 until a stop condition is met, such as the next step (i) leaves the 
domain, (ii) hits the start or end of another ridge, (iii) fails to hnd a maximum or (iv) 
has an FTLE value below a user dehned threshold. 

The fundamental difference between this approach and the previous version is how the FTLE 
held is calculatecP^. To reduce the number of calculations, the previous version sequentially 
calculates the FTLE held at test points, ignoring regions of the held where there are no 
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a 


First Step 


Second Step 




FIG. 3. (a) Initial step of the ridge tracking. The seed point (red) is a local max along the grid Q. 
An initial step direction is set by the blue arrow. The set of test points are black circles with the 
FTLE values as shown in the inset. The corresponding maximum is shown as a red dot in both 
the main and inset plots, (b) Second step in the ridge tracking process. The direction is set by the 
vector from the seed point and the initial step (red dashed arrow). 


FTLE ridges. Because the simultaneous calculation of the FTLE held is more efficient than 
sequential calculation, however, we calculate the FTLE held for the entire domain and use 
interpolation to calculate the ETLE values for the above algorithm. This proves sufficient 
to hnd a reasonable hrst approximation of the ridge, which can then be rehned. 

For systems where the FTLE held has sharp ridges, and thus there are large magnitudes 
of the second-derivative in the direction normal to the ridge, it is more appropriate to select 
an initial step in the direction normal to the gradient. This is due to large gradients normal 
to the ridge for positions not precisely on the ridg^. For subsequent steps, the approach 
does not rely on the gradient calculation, so the method accurately tracks the FTLE ridges. 
Regardless of ridge sharpness, errors in the ridge tracking step are particularly large in the 
initial steps due to not being directly on the ridge. This error is reduced by the rehnement 
approach discussed next. 

Having used the above recipe to initially locate a ridge, a rehnement scheme is then 
employed. As depicted in £gure|^a), the rehnement scheme takes an initial ridge (red), 
calculates the normal at all points, and then places a number of test points at incremental 
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distances 6 either side of the ridge, along the normal direction (black points). FTLE values 
for all these points normal to the ridge are calculated (not interpolated) and the point with 
the maximum FTLE value at each cross ridge location is taken as the revised position of the 
ridge (green circles). For the updated position of the ridge, the normal along the ridge is 
recalculated, a smaller search range either side of the revised ridge is considered, the FTLE 
values for points in this range are calculated, and the ridge position is updated, as illustrated 
in £gure|^b). This process is repeated iteratively for progressively smaller values of 5 until 
the accepted degree of accuracy is reached. 

One of the principal benehts of this method over previous methods is that we hnd that 
an initial calculation of the entire FTLE held is signihcantly faster than running the FTLE 
calculation of a small set of points hundreds of times; for even moderate resolution, the 
interpolated based ridge tracking still produces an acceptable hrst estimate for the ridges of 
the FTLE held, which is the initial guess needed to proceed to the ridge refinement scheme. 
Practical issues to be aware of are: (i) that sharp changes in ridge direction will be difficult 
to resolve because the tangent calculation will be noisy and the transverse search interval of 
consecutive points on the ridge may intersect; and (ii) in situations where there is a complex 
system of closely aligned ridges, the ridge rehnement may jump from one ridge to another 
if the search window overlaps nearby ridges. Typically, these issues can be resolved on a 
case-by-case basis by decreasing the step size in ridge tracking and rehning over a smaller 
window many times. 


B. Results 

To demonstrate the ridge tracking and rehnement approach, we apply these techniques to 
our model system. We used a discretized velocity held with Ax = 2“^^, relative tolerance of 
10“^ for ode45, and cluster size Sa = 10“®. The FTLE held is calculated on a uniform held 
of resolution 6x = 0.01. The ridges found using the ridge tracking method are presented as 
black lines in hgurej^a). These do not smoothly follow the FTLE held due to the underlying 
discretization; this is clearly seen in the ridges near the left and right boundaries and on 
closer inspection of the center ridge. Implementing the rehnement scheme eliminates these 
artihcial features, demonstrated by the green rehned ridges that do smoothly follow the 
FTLE held. 
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FIG. 4. (a) Initial ridge detected via ridge tracking (red) with the points normal to the ridge 
(black). The points closest to the true ridge (blue) have the highest FTLE values and are circled 
(green), (b) The updated position of the estimated ridge using the circled points in (a) is in red. 
Normal points spanning a smaller search range (black) are again calculated. The new maximum 
position is circled in green. 



FIG. 5. (a) Tracked ridges (black), and the corresponding refined ridges (green) of the FTLE field 
(red) for the autonomous system. An inset is included to demonstrate the difference between the 
results, (b) Results of advection for the tracked and refined ridges. The inset highlights how the 
small initial differences grow for the central ridge, (c) FTLE values, ^(s), along the vertical tracked 
and refined ridge as a function of the in-line coordinate, s. 

To demonstrate the importance of refining the FTLE ridge, in hgure |^b) we present 
results for the advection of the unrehned and rehned ridges. The unrehned ridges along 
the left and right curved boundaries are actually advected into the domain, rather than 
remaining on the boundary as they should, and the unrehned ridge across the central region 
of the domain becomes jagged. In contrast, the rehned ridges behave as expected and remain 
smooth. The second beneht of rehnement is presented in hgure |^c). Here the FTLE value 
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is presented for the left ridge as a fnnction of its in-line coordinate, s. There are large 
oscillations in the FTLE value and while the result does appear to be continuous there are 
no features in the system of this scale, so it is clear that these variations are simply a result of 
an insufficiently rehned ridge. Rehnement not only smooths out the results signihcantly, but 
we see that the FTLE ridge values are uniformly higher. Similar results were also obtained 
for the central ridge. 

IV. FTLE RIDGE CLASSIFICATION 

Because the FTLE held only represents the magnitude of the largest eigenvalue of the CG 
tensor, there is ambiguity regarding the deformation associated with an FTLE ridge, which, 
if available, would reveal a great deal about its fate and inhuence. A hrst attempt at FTLE 
classihcation dehned three metrics measuring the normal repulsion, shear, and tangential 
stretching of the ridged^, while these metrics provide some of the desired information, 
the implementation relied on mapping the Hessian of the FTLE forward in time, which is 
numerically challenging. In this section, we therefore further these ideas and modify them 
for application to rehned FTLE ridges. 

A. Methods 

Given a material line, 70 , at initial time, to, for each point along the line, Xq G 70 , 
the tangent and the normal vectors are eg and no, respectively. As presented in hgure 
the material line and its inhnitesimal vectors are mapped forward by the how map. By 
dehnition, the tangent vector remains tangent to the material line, || e*, though 

stretching or contraction of the unit vector is possible. Deformation of the initial normal 
vector, however, is such that it need not remain normal, and so the advected normal can be 
decomposed into two parts: 

VFlno = prit + aCf, (13) 

where rit and e* are the normal and tangent vectors of the advected material line, respec¬ 
tively, p is the normal stretching and a is the Lagrangian sheai^^. 

To classify the inhuence of the FTLE ridges, we consider the relative strengths of normal 
versus tangential growth, and the relative strengths of normal stretching and Lagrangian 


16 


Refining and classifying FTLE ridges 



FIG. 6 . Sketch of the advection of a material line ( 70 ) and the normal (no) and tangent (eg) 
vectors at a point (ico) on the line. 

shear. In this manner, through mapping the vectors Cq and no, we hrst identify if material 
initially along or normal to the ridge is more greatly influenced by the local deformation. 
Formally, we can quantify the magnitudes of the advected normal and tangential vectors 
via. 


ni{xo) = log||VF/^(a;o)no(a;o)||, (14) 

ei{xo) = log||VF/^(a;o)eo(aJo)||, (15) 


where we have added the subscript I to denote the logarithm scaling. Because the loga¬ 
rithm scaling is used, stretching and contraction of the initial unit vectors are represented 
by positive and negative values, respectively. These metrics are related to the previously 
mentioned classihcation schem^®, with the difference being that in this case both measures 
are concerned only with total growth and not their growth relative to the tangent and normal 
directions. In a similar manner, the hyperbolic repulsion and Lagrangian shear are given 
by: 


Pi = log 
(Ti = log 


(rii, VF/^no) 
(eo VF/^no) 


(16) 

(17) 


These metrics do match the previous classihcations, except for a time averaging, the differ¬ 
ence being in how e* and rit are calculated. Note that the positive and negative values of 
Pi correspond to repulsion and attraction in the normal direction, and that the convention 
(i.e. clockwise or counter-clockwise) of shear is lost by the absolute values and only the 
magnitude of shear is measured. The key to calculating these quantities is to accurately 
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identify the ridges so that advection is possible, enabling e* to be reliably calculated by 
applying the flow map gradient to Cq; n* is then found by virtue of it being perpendicular 
to et. 

As we have advocated throughout the paper, in performing any such calculations it is 
important to understand the inherent limitations given that FTLE ridges are the most 
unstable transport features in the flow held. Because the eigenvectors of the how map 
gradient, which we refer to as and ^ 2 , with corresponding eigenvalues A 2 > Ai > 0, 
are orthogonal, they provide a natural basis to decompose the normal and tangential ridge 
vectors and assess the impact of small errors. As demonstrated in appendix when the 
normal and tangential ridge vectors are not closely aligned with ^1 and ^ 2 , small errors 
do not signihcantly impact the classihcation metrics. If the tangent and normal vector are 
approximately aligned with the eigenvector helds, however, small errors become a concern. 
If an FTLE ridge is ehectively a strainline (i.e. everywhere tangent with ^ 1 ) errors do 
not greatly impact n; but can impact ei, though only large amounts of error can cause 
the relative size of these two to be signihcantly altered. Furthermore, for FTLE ridges 
that are close to being strainlines, both pi and ai are sensitive to errors, particularly when 
A 2 S> Ai. Full derivations of the sensitivity of the classihcation metrics that underly this 
summary are presented in appendix This sensitivity analysis is particularly enlightening 
when considering ridges that closely align with the eigenvector helds where the classihcation 
metrics are most sensitive. 


B. Results 


The classihcation metrics (14), (15), (16) and ( p!^ are applied to the analytical model 
system, using the same resolution of discretized data used for the results in hgure The 
classihcation values as a function of a normalized line coordinate, s, are presented in hg- 
ures [^a) and[^b). For this system, the values of Ai and A 2 , and the orientation of the 
ridges relative to the eigenvector held, all lend themselves to reliable classihcation. The 
ridges along the boundary have greater normal than tangential growth along their entire 
length; indeed, only where these boundary ridges are near to the central ridge does their 
tangential growth become positive, and the negative values elsewhere indicate tangential 
contraction. Similarly, there is growth of the normal vector everywhere along the ridge 
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across the center of the domain, and tangential contraction everywhere except near the left 
and right boundary. For all three ridges, the relative amounts of normal repulsion and La- 
grangian shear are approximately equal, with the exception that near the maximum value 
of ni along the boundary ridges, normal repulsion dominates. 

With this information in hand, we now know that as the central ridge is advected there 
will be contraction of material initially distributed along the ridge, stretching of material 
initially normal to the ridge, and that this latter stretching will be evenly distributed between 
shear and normal repulsion. Thus, we expect a circular patch of particles initially placed 
anywhere along the ridge to become a tilted ellipse with its semi-major axis oriented at some 
clearly observable angle to the ridge. This behavior is illustrated by the light green and dark 
blue patches in hgure hguresj^c) and[^d). For the ridges along the boundaries, there is a 
segment where n;, > 0 indicating that there will be stretching both along the tangent and 

normal to the tangent. As demonstrated by the dark green patches in hguresj^c) and[^d), 
we see there is a region of very low Lagrangian shear so that the patch becomes stretched 
normally to the ridge without any signihcant tilting. In contrast, the light blue patch of 
particles released in the vicinity of the intersections between the boundary and central ridges 
experiences both tangential and normal growth; the growth of both components is reflected 
in the non-orthogonal intersection of the patch and ridge and the growth of the patch away 
from the ridge. The only condition not presented in this model is where ei > ni^ which 
would result in an initially circular patch stretching along a ridge without having signihcant 
expansion of material initially normal to the ridge. 


V. APPLICATION TO NINGALOO PENINSULA 

As a further application of our FTLE-based methods, we investigated a data set pro¬ 
duced by a numerical ocean model of a region containing the Ningaloo peninsula in Western 
Australia, which is the location of one of the world’s longest fringing coral reefs and vast 
offshore hydrocarbon reserves. The sea surface velocity helds employed in the study were 
produced by a double nested validated ROMS model with 2000m resolution, and we consider 
a 108-hour time window. A previous LCS analysis of this region investigated the impact of 
incorporating surface wind effects into strainline LCS analysiJ^. We now calculate, advect, 
and classify the FTLE ridges of this system. 
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FIG. 7. (a) Normal (thick lines) and tangential (thin lines) growth along the left (which is sym¬ 
metric to the right) boundary (red liness) and center (black lines) FTLE ridges as a function of a 
normalized inline coordinate, s. (b) Normal repulsion (thick) and Lagrangian shear (thin) along 
the FTLE ridges. Advection of FTLE ridges and passive tracer patches from time (c) t = 0 to (d) 
t = 2. 


Before analyzing the FTLE field, we performed a sequence of convergence studies to 
confirm reliable results. First the relative tolerance of the ODE solver was varied to ensure 
that the flow map was calculated accurately, and the relative tolerance of 10“^ was sufficient 
for converged trajectory calculations. Next, cluster sizes were varied to hnd converged values 
of the FTLE held. Cluster sizes of Im, 10“^m, and 10“‘^m were tested and the results for 
all three cluster sizes were nearly identical, so a cluster size of Im was used. While this 
may seem small given a grid resolution of 2000m, it should be noted that the ratio of grid 
resolution to cluster size is approximately the same as that used for the discretized studies 
of our analytical model. As a hnal test, we subsample the domain to ensure that the results 
are not signihcantly affected. When decreasing the domain resolution to 4000m, we found 
that the dominant features of our analysis were nearly identical, with only slight (0(5%)) 
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FIG. 8. (a) Forward-time (light and dark red) and the backward FTLE ridges (light blue) along 
with passive tracer patches (light and dark green) at times to =0:00 21 December 2009 and ti = 
12:00 25 December 2009 overlaying a MODIS satellite image of the Ningaloo peninsula. Inset 
of Australia shows the location of the region. Positive and negative shearlines (black and gray), 
strainlines (yellow), and (b) forward-time FTLE ridges or (c) backward-time FTLE ridges. 

differences in absolute value of the FTLE maxima. Having determined the parameters 
needed for a robust FTLE calculation, the forwards- and backwards-time FTLE helds on a 
200 m grid were calculated in order to perform the initial ridge extractions, followed by ridge 
rehnement and classihcation. 

Two forward-time ridges that are in close proximity and a single backward time ridge 
were particularly strong features of the FTLE held, and so we focus our attention on these. 
The positions of all three ridges at the beginning and end of the time window, which required 
advection of the rehned ridges, are presented in hgurej^a). Particle patches have been added 
to help demonstrate the local deformations near the ridges. A Lagrangian hyperbolic point 
lies at the intersection of one of the forward-time FTLE ridges and the backward-time FTLE 
ridge. The other forward time FTLE ridge, however, was not associated with any backward 
time ridge. Figures |^b) and (c) plot strainlines and shearlines in the vicinity of these FTLE 
ridges, and we see that the two FTLE ridges that form the Lagrangian hyperbolic point are 
nicely captured by strainlines; the isolated forward-time FTLE ridge, however, is not well 
captured by a single strainline or shearline. 

The classihcation scheme is implemented on these three ridges, and the values of ei and 
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ni(s) 


e/(s) 




FIG. 9. (a) n; (bold) and e; (thin) as a function of s for the isolated forward-time ridge (light red) 
and intersecting forward-time ridge (dark red), (b) n; and ei along the backwards-time ridge. 


rii are presented as a function of a dimensionless inline coordinate s in figure [^a) and [^b) 
(s = 0 corresponds to the right most end of the ridge). For all three ridges there is growth of 
the normal vector along the entire length of the ridge, indicating that the ridge is influencing 
fluid on either side. The persistent negative values of e; along the backward-time ridge means 
that this ridge will stretch in tangential length from time to to ti (contraction in backwards 
time). The ei value of the forward-time ridge that intersects the backwards time ridge is 
also persistently negative indicating that it will contract towards the backward-time ridge. 
Interestingly, the third FTLE ridge has regions where < 0 and ei > 0, indicating there are 
alternating segments of expansion and contraction. Additionally, we note that the values of 
ei along the two intersecting ridges appears to be noisy while the isolated forward-time ridge 
has a smoothly varying ei. This is indicative of the fact that the two intersecting ridges are 
well represented by strainlines. 

Finally, the values of pi and ai are investigated along each of the three ridges with the 
results presented in hgure It should be noted that in all three cases the value of ni is 
smoothly varying along the entire length of the line (presented in black for reference). The 
intersecting forward and backward-time ridges, presented in [l^a) and [T^b), respectively 
demonstrate large variation in the values of pi and ai along the line. This error is a mani¬ 
festation of the sensitivity of these classihcation metrics when a ridge closely aligns with the 
CG eigenvector held (see Appendix B for further details). The smoothly varying pi and ai 
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FIG. 10. Pi (green), ai (blue), and ni (black) as a function of the scaled inline coordinate, s for 
the (a) intersecting forward-time ridge, (b) intersecting backward-time ridge, and the (c) isolated 
forward-time ridge. 

of the isolated ridge are presented in figure [T^c). While there are a couple isolated points 
of small Lagrangian shear, most of the ridge features larger values of Lagrangian shear than 
hyperbolic repulsion. This manifests itself in the tilting of the particle patch relative to the 
ridge in hgure These results are more reliable than the intersecting ridge classihcations 
because in this instance, the ridges are not closely aligned with the eigenvector helds placing 
it in a regime where the metrics are robust to small errors. 


VI. CONCLUSIONS 

While the FLTE method has its known limitations, it nevertheless identihes a particular 
type of LCS, namely the most active material lines in a fluid flow, and is a useful tool 
for preliminary analysis of flow transport. Of particular importance, therefore, is that the 
limitations of what can be interpreted from the FTLE held are clear and that proper care 
is taken to ensure that the FTLE held and its ridges are accurately determined. 

In this paper, we used both the hnite-diherence and advected-gradient approaches for 
calculating the FTLE held, and found that any benehts of the advected-gradient method 
are lost when working with discretized data. While it is reassuring that the hnite-diherence 
method is found to be robust, convergence of the how map by varying the relative tolerance 
and of the hnite-diherences by varying cluster size are necessary. The normal-maximum 
ridge dehnition was introduced as a combination of the strengths of both the height and 
watershed ridges. Recognizing the ridge detection scheme by Lipinski and Mohsenil^^, minor 
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modifications were made to take advantage of the parallelization of Matlab; a practical 
scheme was introduced to enable rapid rehnement of an FTLE ridge. 

The rehned FTLE ridges mark material lines with large associated deformation, but 
classihcation and advection is furthermore necessary to understand how these material lines 
are influencing flow transport. Measuring the growth of the normal and tangent vectors 
determines if an FTLE ridge is primarily influencing material along its length or adjacent 
to it. In cases where normal growth dominates, it is then reasonable to compare the normal 
repulsion and Lagrangian shear. We performed some analysis to determine when these 
different types of deformation are sensitive to errors in the FTLE ridge. As a reasonable 
test, these techniques were applied to both an analytic model and an ocean surface data set. 
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Appendix A: Analytic model 

The following is a basic analysis of the autonomous system in order to get the FTLE held. 
In order to calculate the FTLE held, it is necessary to obtain the how map and its gradient 
in terms of the transformed coordinates. First we can solve equations ([^ by separation of 
variables yielding: 


Xi{t,Ai) = sign(Ai) 
A2(f, As) = sign(A2) 


1 

1 



(Al) 


where (Ai, As) are the Lagrangian coordinates, i.e. initial positions, of the trajectories. Then 


to get the how map of the transformed system, we use substitution of the (Al) into (10). 
Because it is convenient to present the how map in terms of its initial conditions in the 
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transformed coordinates, we note that 


Ai = ai cos(ro) + 02 sin(ro), 
A 2 = 02 cos(ro) - Oi sin(ro), 


(A2) 


where (oi, 02) are the Lagrangian coordiantes in the transformed system and tq = 

-y/oi^ + 02^. 

it 


Next to get the flow map gradient, V-F/^ = dxi/daj, the chain rule is applied to (10) 


Taking the derivative of equation (10) in terms of the Lagrangian coordinate a* gives: 
dxi 


d 

= — [Xi cos(r) - X2 sin(r)] , 


dai da. 

'dXi 


-X — 
dai dai 


cos r — 


dai ^ dai 


sin(r). 


f)T f) 

^ = — [X2 cos(r) + Xi sin(r)] , 


dai dai 

'dX2 


dai ^ dai 


cos(r) + 


dXi 

dai 




sm r . 


(A3) 


(A4) 


Equations (A3) and (A4) indicate that we need to start deriving the terms dXi/daj and 


dr/dai- From solving equation (10) with the given initial conditions, it can be shown that 


dX, 

dAi 


0 


0 

which can be utilized in the chain rule expansion to show that: 

dX, ^ dAk dXj 
daj daj dA^ 

If we reconhgure the coordinate transform it can be shown that 

Ai = ai cos(ro) — 02 sin(ro), 

A 2 = 02 cos(ro) + Oi sin(ro), 

where tq = \/a^^ + 02^. This leads to 


(AS) 


(A6) 


dai 


A 

ro 
, 2 


1 - ^ ) cos(r) - ^ sin(ro) - (1 + ^ ) sin(r) - ^ cos(ro) 


ro 


1 - ^ ) sin(r) + ^ cos(ro) (1 + ^ ) cos(r) - ^ sin(ro) 


ro 


ro 


ro 


(A7) 


which when substituted along with equation (AS) give us the necessary terms to calculate 


equation (A6). 
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Calculating the remaining derivatives of equations (A3) and (A4) is via 


dr 

dtti 


dXj dr 
dtti dXi 


d^jq 

dtti r 


(AS) 


where 


dr 


r 


With the necessary derivatives to calculate the terms of the flow map gradient, we can 
now calculate the terms of the CG tensor: 


Clixo) = 


Gn Gi2 


Gi2 C 22 



dxi ^ I dx2 ^ dxi dx\ i dx2 dx2 

da\ ' dai dai da2 ~ da\ da2 

dxi dxi I dx2 8x2 dxi ^ i 8x2 ^ 

9ai da2 ' dai da2 da2 ' da2 


(A9) 


Finally, with the generic terms of the CG tensor, the largest eigenvalue and eigenvector can 
be shown to be: 


Ao — 


Cii + C22 + 


(AlO) 


Appendix B: Sensitivity analysis 

Central to the following arguments is the singular value decomposition of the flow map 
gradient. We dehne the smaller singular value as Ai and it has the corresponding right- 
singular vector, ^1, which maps forward to the left-singular vector, ui such that 

= Ai^i, 

= ^2^2- 

Throughout the derivation we utilize the conventions $,i x ^2 = x U 2 = 1. 

A point, Xq, along the FTLE ridge, 70, is allowed to have the general orientation: 

Co = a/1 - + b$,2 

^0 = + a /1 - b ‘^^ 2 - 

where |6| < 1 and the tangent and normal vectors have unit length. In the unique case 
where the FTLE ridge is a strainline (stretchline), b = 0 {b = 1). For the following analysis, 
a perturbation e 1 is added to the tangent vector such that the ^2 coefficient becomes 
b + e, and so the other components of this vector and of the normal vector are modihed to 
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that both remain of unit length to leading order. This leads to the perturbed tangent and 
normal vectors: 

eo = (v/l - "'2(1-^2)372 + + (^ + e)6, 

< = -{b + e)6 + (v^l - - "'2(1 -^ 2 ) 3/2 + ^2- 

where ' denotes the perturbed vector. Mapping these forward under the action of the flow 
map gradient gives: 

- e-j=== - e^2(i _ ^2)3/2 + + (^ + e)A2W2, 

= -{b + e)\iUi + - 6^ - e-j== - <^^2(1 - 52)3/2 + -^2^i2, 

and it can be shown that: 

e[ = A2 [ 5 " + {b + ef{l - 6 ^) + 5 ^ 0 {e^)\ , 

n'l = A2 [1 - (6 + efil - 5 ') + 0(e=')] . 


The advected tangent and normal are thus: 


L 


y/l^ 


vr^ 


^t= A7 


n L 


-(6 + e)A 2 Wi + ( - 


2(1 _ 52 ) 3/2 

b 


+ 0(e ) j Aiiti + (6 + e)X2U2 

1 


— e 


Vl^ 2(1 - 62)3/2 


+ 0 (e'*) A 1 W 2 


The normally hyperbolic repulsion and Lagrangian shear can be determined via: 


p’ = {n’„ VFl,n’„) = ^ [1 + 0(£=>)] , 


a' = {e'„ VF‘ni> = 




bVl^ + e 


262 


36 - 26^ 


2(1 - 62)3/2 


+ 0{e^) 


To gain insight into the stability of these derivations to perturbation we need to expand 
the terms and look at the coefficient multiplying e. Of particular interest are three cases. The 
hrst case is strainlines (6 = 0); these correspond to material lines for which normal repulsion 
is maximized and thus are tangent to the ^i-£eld. We also consider the stretchlines (6 = 1), 
for which stretching is principally in a tangential direction. Finally, we consider cases where 
Ai A 2 corresponding to 6 1; this is a common scenario for systems with large amounts 

of stretching. 
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The first classification metric we look at is the growth of the tangent vector, e;. The 
growth of the perturbed tangent vector is obtained from; 


e; = A2[<52 + {b + e)2(l - 52) + <52Q(e3)]V2 


= \2[S^ + b\l-S^)Y/^ + ^ 


6(1 - 52) 


v/52 + 62(1 _ 52) 


+ 0(6^) 


The coefficient for the perturbation remains small for all cases we consider. It goes to 1 — 6^ 
as 6 —)■ 1; therefore, c; for a stretchline is not signihcantly altered when the tangent vector 
orientation is perturbed. For 6 1 , the perturbation coefficient goes to 1 , so when the 

vector is not closely aligned to either eigenvector held, perturbations will not result in large 
changes in ei. The perturbation coefficient becomes zero, however, when 6 = 0 , and so we 
expand to consider the second order perturbation. It can be shown that when 6 = 0 , the 
second order expansion becomes: 



This shows that when 6 1 and 6 = 0 , the coefficient for the perturbation diverges unless 

€ <^ 6 . This is not surprising since, for a strainline, the tangent is aligned with the direction 
of smallest stretching, ^1, and the normal is aligned with the direction of largest stretching, 
^2- Any perturbation to these vectors adds a contribution from ^2 to the tangent vector; 
this component grows signihcantly causing the overall length of the advected tangent vector 
to be much larger than the unperturbed state. 

Next we consider the growth of a perturbed normal vector: 


n' = A2[1 - (6 + e)2(l - 6^) + 

For the stretchline case (6 = 1 ), the perturbation coefficient reduces to 1/5 — 5 . When 6 1 , 

this is very large indicating that ni for stretchlines is sensitive to small perturbations; this 
can be rationalized in much the same way that sensitivity of e; for strainlines was. The 
perturbation term again reduces to zero for the strainline case (6 = 0), so the expansion to 
second order is again performed: 


n'l = X2 


1 + h=(l - 
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In this case, small perturbations do not significantly alter the ni value because while the 
absolute error may be large, relative to the unperturbed value it is not. For values of b not 
near the extremes, the perturbation coefficient is of order one meaning that the values are 
not sensitive to perturbation. 

The sensitivity of the normal hyperbolic growth depends on the expansion of l/e'p 


A1A2 



A1A2 r 

A2[<52 + 62(1 _ 52)]1/2 [ ^52 + ^2(1 _ ^2)_ 

For (5 1 and 6 7^ 0 the perturbation coefficient remains small. In the case where 6 = 0 , 

we have shown that is sensitive to perturbation so it is not surprising that p'l will also be 
sensitive to errors. Through substitution, it can be shown that for 6 = 0 : 

Pi = A2 (1 + - 1 )^ . 

Unless e 6 this term is also sensitive to perturbations. 

Finally, the analysis of the Lagrangian shear produces similar results to the hyperbolic 


repulsion. 


Ai-A? 


6 Vr^ + e- 


2 36 - 26^ 3 

^ 2(1 - 62 ) 3/2 + ^ 


+ 0(e^) 


^ (A2 _ Xl)bVT^ r / 1 - 26 ^ _ 6(1 - 6^) A 2; 

A 2[62 + 62 ( 1 - 62 )]i /2 [ + \b{l-b^) 6 ^ + b^{l- 6 ^)) ^ [ 

For 6 = 1 , there is a singularity in the perturbation coefficient indicating that perturbations 

for a stretchline will result in large changes in (t(. Again for the strainline, it is necessary 

to expand to the second order to study the stability of perturbations. This results in the 


expansion: 



It makes sense that this term is only insensitive to noise when e 6, because it relies on e'l 


being insensitive to noise. 
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